home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-19 / iritsm3s.zip / MSHPLANR.C < prev    next >
C/C++ Source or Header  |  1991-11-30  |  10KB  |  306 lines
  1. /******************************************************************************
  2. * MshPlanr.c - Test colinearity of control meshes/polygons.              *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Sep. 91.                          *
  5. ******************************************************************************/
  6.  
  7. #include "cagd_loc.h"
  8.  
  9. /*****************************************************************************
  10. * Computes a distance between two control points at given indices Index?.    *
  11. * Only E2 or E3 point types are supported.                     *
  12. *****************************************************************************/
  13. CagdRType CagdDistanceTwoCtlPts(CagdRType **Points, int Index1, int Index2,
  14.                             CagdPointType PType)
  15. {
  16.     switch (PType) {
  17.     case CAGD_PT_E2_TYPE:
  18.         return sqrt(SQR(Points[1][Index1] - Points[1][Index2]) +
  19.             SQR(Points[2][Index1] - Points[2][Index2]));
  20.     case CAGD_PT_E3_TYPE:
  21.         return sqrt(SQR(Points[1][Index1] - Points[1][Index2]) +
  22.             SQR(Points[2][Index1] - Points[2][Index2]) +
  23.             SQR(Points[3][Index1] - Points[3][Index2]));
  24.     default:
  25.         FATAL_ERROR(CAGD_ERR_UNSUPPORT_PT);
  26.         return 0.0;
  27.     }
  28. }
  29.  
  30. /*****************************************************************************
  31. * Fits a plane into the four points from Points indices Index?. Points may   *
  32. * be either E2 or E3 only.                             *
  33. *   Returns 0.0 if failed to fit a plane, otherwise a measure on the size of *
  34. * the mesh data (distance between points) isreturned.                 *
  35. *****************************************************************************/
  36. CagdRType CagdFitPlaneThruCtlPts(CagdPlaneStruct *Plane, CagdPointType PType,
  37.                      CagdRType **Points,
  38.                      int Index1, int Index2, int Index3, int Index4)
  39. {
  40.     int i, j, Indices[4];
  41.     CagdRType SizeMeasure,
  42.     MaxDist = 0.0;
  43.     CagdVecStruct V1, V2, V3;
  44.  
  45.     Indices[0] = Index1;
  46.     Indices[1] = Index2;
  47.     Indices[2] = Index3;
  48.     Indices[3] = Index4;
  49.  
  50.     /* Find out the pair of vertices most seperated: */
  51.     for (i = 0; i < 4; i++)
  52.     for (j = i + 1; j < 4; j++) {
  53.         CagdRType
  54.         Dist = CagdDistanceTwoCtlPts(Points, Indices[i], Indices[j],
  55.                                      PType);
  56.         if (Dist > MaxDist) {
  57.         MaxDist = Dist;
  58.         Index1 = i;
  59.         Index2 = j;
  60.         }
  61.     }
  62.  
  63.     if (MaxDist < EPSILON) return 0.0;
  64.     SizeMeasure = MaxDist;
  65.     MaxDist = 0.0;
  66.  
  67.     /* Find a third most seperated than the selected two. */
  68.     for (i = 0; i < 4; i++)
  69.     if (i != Index1 && i != Index2) {
  70.         CagdRType
  71.         Dist1 = CagdDistanceTwoCtlPts(Points, Indices[Index1],
  72.                           Indices[j], PType),
  73.         Dist2 = CagdDistanceTwoCtlPts(Points, Indices[Index2],
  74.                           Indices[j], PType),
  75.         Dist = MIN(Dist1, Dist2);
  76.  
  77.         if (Dist > MaxDist) {
  78.         MaxDist = Dist;
  79.         Index3 = i;
  80.         }
  81.     }
  82.     if (MaxDist < EPSILON) return 0.0;
  83.  
  84.     /* Now we have the 3 most seperated vertices to fit a plane thru. */
  85.  
  86.     switch (PType) {
  87.     case CAGD_PT_E2_TYPE:
  88.         /* This is trivial since the plane must be Z = 0. */
  89.         Plane -> Plane[0] = 0.0;
  90.         Plane -> Plane[1] = 0.0;
  91.         Plane -> Plane[2] = 1.0;
  92.         Plane -> Plane[3] = 0.0;
  93.         break;
  94.     case CAGD_PT_E3_TYPE:
  95.         /* Compute normal as a cross product of two vectors thru pts. */
  96.         for (i = 0; i < 3; i++) {
  97.         j = i + 1;
  98.         V1.Vec[i] = Points[j][Index2] - Points[j][Index1];
  99.         V2.Vec[i] = Points[j][Index3] - Points[j][Index2];
  100.         }
  101.         CROSS_PROD(V3.Vec, V1.Vec, V2.Vec);
  102.         CAGD_NORMALIZE_VECTOR(V3);
  103.         for (i = 0; i < 3; i++)
  104.         Plane -> Plane[i] = V3.Vec[i];
  105.  
  106.         Plane -> Plane[3] = (-(V3.Vec[0] * Points[1][Index1] +
  107.                    V3.Vec[1] * Points[2][Index1] +
  108.                    V3.Vec[2] * Points[3][Index1]));
  109.         break;
  110.     default:
  111.         FATAL_ERROR(CAGD_ERR_UNSUPPORT_PT);
  112.         break;
  113.     }
  114.  
  115.     return SizeMeasure;
  116. }
  117.  
  118. /*****************************************************************************
  119. * Computes and returns distance between point Index and given plane which is *
  120. * assumed to be normalized, so that the A B C D plane normal has unit len..  *
  121. * Also assumed the Points are non rational with MaxDim dimension.         *
  122. *****************************************************************************/
  123. CagdRType CagdDistPtPlane(CagdPlaneStruct *Plane, CagdRType **Points,
  124.                               int Index, int MaxDim)
  125. {
  126.     int i;
  127.     CagdRType
  128.     R = Plane -> Plane[3];
  129.  
  130.     for (i = 0; i < MaxDim; i++)
  131.     R += Plane -> Plane[i] * Points[i + 1][Index];
  132.  
  133.     return fabs(R);
  134. }
  135.  
  136. /*****************************************************************************
  137. * Computes and returns distance between point Index and the line from first  *
  138. * point in direction LineDir (usually the line direction to last the point). *
  139. * LineDir is assumed to be unit length normalized vector.             *
  140. *****************************************************************************/
  141. static CagdRType CagdDistPtLine(CagdVecStruct *LineDir, CagdRType **Points,
  142.                         int Index, int MaxDim)
  143. {
  144.     int i;
  145.     CagdRType R;
  146.     CagdVecStruct V1, V2;
  147.  
  148.     for (i = 0; i < MaxDim; i++)
  149.     V1.Vec[i] = Points[i+1][Index] - Points[i+1][0];
  150.     if (MaxDim < 3)
  151.     V1.Vec[2] = 0.0;
  152.  
  153.     /* Let V1 be the vector from point zero to point index. Then the distance */
  154.     /* from point Index the the line is: | (LineDir . V1) LineDir - V1 |.     */
  155.     V2 = *LineDir;
  156.     R = DOT_PROD(V1.Vec, V2.Vec);
  157.     CAGD_MULT_VECTOR(V2, R);
  158.     CAGD_SUB_VECTOR(V2, V1);
  159.  
  160.     return CAGD_LEN_VECTOR(V2);
  161. }
  162.  
  163. /*****************************************************************************
  164. * Test polygon colinearity by testing the distance of interior control       *
  165. * points from the line connecting the two polygon end points.             *
  166. *   Returns a relative ratio of deviation from line relative to its length.  *
  167. *   Zero means all points are colinear.                         *
  168. *   If two end points are same ( no line can be fit ) INFINITY is returned.  *
  169. *****************************************************************************/
  170. CagdRType CagdEstimateCrvColinearity(CagdCrvStruct *Crv)
  171. {
  172.     int i,
  173.     MaxDim = 3,
  174.     Length = Crv -> Length,
  175.     Length1 = Length - 1;
  176.     CagdRType LineLen,
  177.     MaxDist = 0.0,
  178.     **Points = Crv -> Points;
  179.     CagdCrvStruct
  180.     *CoercedCrv = NULL;
  181.     CagdPointType
  182.     PType = Crv ->PType;
  183.     CagdVecStruct LineDir;
  184.  
  185.     switch (PType) {          /* Convert the rational cases to non rational. */
  186.     case CAGD_PT_P2_TYPE:
  187.         CoercedCrv = CagdCoerceCrvTo(Crv, CAGD_PT_E2_TYPE);
  188.         Points = CoercedCrv -> Points;
  189.         PType = CoercedCrv -> PType;
  190.         break;
  191.     case CAGD_PT_P3_TYPE:
  192.         CoercedCrv = CagdCoerceCrvTo(Crv, CAGD_PT_E3_TYPE);
  193.         Points = CoercedCrv -> Points;
  194.         PType = CoercedCrv -> PType;
  195.         break;
  196.     }
  197.  
  198.     switch (PType) {
  199.     case CAGD_PT_E2_TYPE:
  200.         LineDir.Vec[0] = Points[1][Length1] - Points[1][0];
  201.         LineDir.Vec[1] = Points[2][Length1] - Points[2][0];
  202.         LineDir.Vec[2] = 0.0;
  203.         MaxDim = 2;
  204.         break;
  205.     case CAGD_PT_E3_TYPE:
  206.         LineDir.Vec[0] = Points[1][Length1] - Points[1][0];
  207.         LineDir.Vec[1] = Points[2][Length1] - Points[2][0];
  208.         LineDir.Vec[2] = Points[3][Length1] - Points[3][0];
  209.         break;
  210.     default:
  211.         FATAL_ERROR(CAGD_ERR_UNSUPPORT_PT);
  212.         break;
  213.     }
  214.  
  215.     LineLen = CAGD_LEN_VECTOR(LineDir);
  216.     if (LineLen < EPSILON) {
  217.     if (CoercedCrv != NULL)
  218.         CagdCrvFree(CoercedCrv);
  219.     return INFINITY;
  220.     }
  221.  
  222.     CAGD_DIV_VECTOR(LineDir, LineLen);
  223.  
  224.     for (i = 1; i < Length1; i++) {
  225.     CagdRType
  226.         Dist = CagdDistPtLine(&LineDir, Points, i, MaxDim);
  227.  
  228.     if (Dist > MaxDist)
  229.         MaxDist = Dist;
  230.     }
  231.  
  232.     if (CoercedCrv != NULL)
  233.     CagdCrvFree(CoercedCrv);
  234.  
  235.     return MaxDist / LineLen;
  236. }
  237.  
  238. /*****************************************************************************
  239. * Test mesh colinearity by testing the distance of interior points from the  *
  240. * plane thru 3 corner points.                             *
  241. *   Returns a relative ratio of deviation from plane relative to its size.   *
  242. *   Zero means all points are coplanar.                         *
  243. *   If end points are same ( no plane can be fit ) INFINITY is returned.     *
  244. *****************************************************************************/
  245. CagdRType CagdEstimateSrfPlanarity(CagdSrfStruct *Srf)
  246. {
  247.     int i,
  248.     ULength = Srf -> ULength,
  249.     ULength1 = ULength - 1,
  250.     VLength = Srf -> VLength,
  251.     VLength1 = VLength - 1;
  252.     CagdRType
  253.     PlaneSize = 0.0,
  254.     MaxDist = 0.0,
  255.     **Points = Srf -> Points;
  256.     CagdSrfStruct
  257.     *CoercedSrf = NULL;
  258.     CagdPointType
  259.     PType = Srf ->PType;
  260.     CagdPlaneStruct Plane;
  261.  
  262.     switch (PType) {          /* Convert the rational cases to non rational. */
  263.     case CAGD_PT_P2_TYPE:
  264.     case CAGD_PT_E2_TYPE:
  265.         /* Trivial case - it is planar surface. */
  266.         return 0.0;
  267.     case CAGD_PT_P3_TYPE:
  268.         CoercedSrf = CagdCoerceSrfTo(Srf, CAGD_PT_E3_TYPE);
  269.         Points = CoercedSrf -> Points;
  270.         PType = CoercedSrf -> PType;
  271.         break;
  272.     }
  273.  
  274.     switch (PType) {
  275.     case CAGD_PT_E3_TYPE:
  276.         PlaneSize = CagdFitPlaneThruCtlPts(&Plane, PType, Points,
  277.                     CAGD_MESH_UV(Srf, 0,        0),
  278.                     CAGD_MESH_UV(Srf, 0,        VLength1),
  279.                     CAGD_MESH_UV(Srf, ULength1, 0),
  280.                     CAGD_MESH_UV(Srf, ULength1, VLength1));
  281.         break;
  282.     default:
  283.         FATAL_ERROR(CAGD_ERR_UNSUPPORT_PT);
  284.         break;
  285.     }
  286.  
  287.     if (PlaneSize < EPSILON) {
  288.     if (CoercedSrf != NULL)
  289.         CagdSrfFree(CoercedSrf);
  290.     return INFINITY;
  291.     }
  292.  
  293.     for (i = ULength *VLength; i > 0; i--) {
  294.     CagdRType
  295.         Dist = CagdDistPtPlane(&Plane, Points, i, 3);
  296.  
  297.     if (Dist > MaxDist)
  298.         MaxDist = Dist;
  299.     }
  300.  
  301.     if (CoercedSrf != NULL)
  302.     CagdSrfFree(CoercedSrf);
  303.  
  304.     return MaxDist / PlaneSize;
  305. }
  306.